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ABSTRACT 

We test the hypothesis that metal-poor globular clusters form within disk galaxies at redshifts z > 3. Nu- 
merical simulations demonstrate that giant gas clouds, which are cold and dense enough to produce massive 
star clusters, assemble naturally in hierarchical models of galaxy formation at high redshift. Do model clusters 
evolve into observed globular clusters or are they disrupted before present as a result of the dynamical evo- 
lution? To address this question, we calculate the orbits of model clusters in the time-variable gravitational 
potential of a Milky Way-sized galaxy, using the outputs of a cosmological A^-body simulation. We find that 
at present the orbits are isotropic in the inner 50 kpc of the Galaxy and preferentially radial at larger distances. 
All clusters located outside 10 kpc from the center formed in satellite galaxies, some of which are now tidally 
disrupted and some of which survive as dwarf galaxies. The spatial distribution of model clusters is spheroidal 
and the fit to the density profile has a power-law slope 7 2.7, somewhat shallower than but consistent with 
observations of metal-poor clusters in the Galaxy. The combination of two-body relaxation, tidal shocks, and 
stellar evolution drives the evolution of the cluster mass function from an initial power law to a peaked distri- 
bution, in agreement with observations. However, not all initial conditions and not all evolution scenarios are 
consistent with the observed mass function of the Galactic globular clusters. The successful models require the 
average cluster density, M/R 3 h , to be constant initially for clusters of all mass and to remain constant with time. 
Synchronous formation of all clusters at a single epoch (z = 4) and continuous formation over a span of 1 .6 Gyr 
(between z = 9 and z = 3) are both consistent with the data. For both formation scenarios, we provide online 
catalogs of the main physical properties of model clusters. 

Subject headings: galaxies: formation — galaxies: kinematics and dynamics — galaxies: star clusters — 
globular clusters: general 
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1. INTRODUCTION 

Observations in the past decade with the Hubble Space 
Telescope and major ground-based telescopes have revolu- 
tionized the field of star cluster research. The classical di- 
vision into old, massive globular clusters and young, dif- 
fuse open clusters has been challenged by the discoveries of 
young massive star clusters in nearby merging, interacting, 
and even normal star forming galaxies (e.g., iHoltzma neLalJ 
1 19921 IWhitmore et all fi993l IWhitmore & Schweizeilll995t 
b'Connell et all 1 9951 IZepfet alJll999HLarseiJl2002Tr These 
young clusters have the range of masses and sizes similar to 
Galactic globulars and are believed to evolve eventually into 
globular clusters (e.g., [H o & Filinpenko 1996; Meng elet alJ 
120021 Ide Griis et alJl2004l) . However, the mass functions of 
old and young clusters are very different. While the mass 
distribution of globular clusters has a well-defined peak at 
w2x 10 5 M Q , the mass function of young clusters is con- 
sistent wifhasir^epowerlaw between 10 4 and 10 7 M f 7i 



4?. 



(Zhan g&Falll999tldeGriis et alJ2003UAnders et alJ20O 

Did old globular clusters form by a similar physical process 
to the young clusters but gradually evolve into the observed 
distribution ? 

This question is vital to understandi ng early galaxy for- 
mation and is still a matter of debate (Fall & Zhand l2~001t 
Vespe rini et al J 120031) . Globular clusters contain relic infor- 
mation about star formation processes in galaxies at high 
redshift and have been often used to constrain cosmologi- 
cal theory. Uncovering the formation of primeval galaxies 
is a major driving force behind the development of future 
major ground and space observatories. In contrast to the 
formidable difficulties of detecting the birthplaces of glob- 



ular clusters at high redshift, local and interacting galax- 
ies (M33, M51, M82, the Antennae) offer plenty of detail 
of the formation of young star clusters in giant molecular 
cloud s JWilson et all2003tlEng ar giola et all2003UKeto et alJ 
2005). Using ultrahigh-resolution cosmological simula- 
tion, [Kravtsov & Gnedin (2005) identified similar molecular 
clouds within their simulated high-redshift galaxies and pro- 
posed that they harbor massive star clusters. In this paper, 
we investigate whether the dynamically-evolved distribution 
of these model clusters can be reconciled with the current 
p roperties o f globular c lusters . 

Kravtsov & Gnedin (2005) find that giant clouds assem- 
ble during gas-rich mergers of progenitor galaxies, when the 
available gas forms a thin, cold, self-gravitating disk. The 
disk develops strong spiral arms, which further fragment into 
separate molecular clouds located along the arms as beads 
on a string (see their Fig. 1). As a result, massive clus- 
ters form in relatively massive galaxies, with the total mass 
Mhaio > 10 9 Mq, beginning at redshift z ~ 10. The mass and 
density of the molecular clouds increase with cosmic time, 
but the rate of galaxy mergers declines steadily. Therefore, the 
cluster formation efficiency peaks at a certain extended epoch, 
around z ~ 4, when the Universe is only 1 .5 Gyr old. The mass 
function of model clusters is consistent with a power law 



dN I dM = NqM~° , 



(1) 



where a = 2.0 ±0.1, similar to the observations of nearby 
young star clusters. According t o the prescription fo r metal 
enrichment by supernovae in the Kravtsov & Gnedinl ( 120051) 
simulation, model clusters have iron abundances [Fe/H] < — 1 
at z > 3. They correspond to the metal-poor (dominant) sub- 
population of the Galactic globular clusters, which we will 
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use for the comparison. 

We adopt this model to set up the initial positions, veloci- 
ties, and masses for our globular clusters. We then calculate 
cluster orbits using the results of a separate collisionless N- 
body simulation described in § 121 This is necessary because 
the original gasdynamics simulation was stopped at z ~ 3, due 
to limited computational resources. By using the A^-body sim- 
ulation of a similar galactic system, but complete to z = 0, 
we are able to follow the full dynamical evolution of globular 
clusters until the present epoch. We use the evolving prop- 
erties of all progenitor halos, from the outputs with a time 
resolution of ~ 10 8 yr, to derive the gravitational potential in 
the whole computational volume at all epochs. We convert a 
fraction of the dark matter mass into the analytical flattened 
disks, in order to model the effect of baryon cooling and star 
formation on the galactic potential. We calculate the orbits 
of globular clusters in this potential from the time when their 
host galaxies accrete onto the main (most massive) galaxy. 
Using these orbits, we calculate the dynamical evolution of 
model clusters, including the effects of stellar mass loss, two- 
body relaxation, tidal truncation, and tidal shocks. Since the 
efficiency of each of these processes depends on both cluster 
mass, M, and half-mass radius, /?/,, we consider several evo- 
lutionary dependencies Ri,(M). We also consider two possible 
formation scenarios, one with all clusters forming in a short 
interval of time around redshift z = 4, and the other with a 
continuous formation of clusters between z = 9 and z = 3. 

Note that we study only the collisionless evolution of star 
clusters, using the analytical formalism developed in previ- 
ous studies and described in § [3] We do not include the 
recently suggested "infant mortality" effect, which may dis- 
solve a large number of unbound clusters within ^100 Myr 
of their formation in dependent of cluster mass (e.g.. lFall et alJ 
120051 iMengel et aTll2005l) . Such universal reduction in num- 
bers is therefore absorbed in our normalization of the initial 
GCMF. The models calculated in this paper apply only to the 
initially bound clusters. 

2. ORBITS OF GLOBULAR CLUSTERS IN A 
HIERARCHICALLY-FORMING GALAXY 

2.1. Cosmological Simulation 

We set the initial conditions and calculate the orbits of 
model clusters using a cosmological A^-body simulation per- 
formed with the Adaptive Refinem ent Tree code (ART, 
iKravtsov et aljfl997l lKravtsovlll999l) . This simulation fol- 
lows the hierarchical formation of a Milky Way-sized halo 
in a comoving box of 25/T 1 Mpc in the concordance ACDM 
cosmology: (l m = 0.3, J7a = 0.7, <r% = 0.9, Ho = 70 km s" 1 
Mpc" 1 . The simulation begins with a uniform 256 3 grid cov- 
ering the entire computational volume. High mass resolution 
is achieved around collapsing structures by recursive refine- 
ment of this regions using an adaptive refinement algorithm. 
The grid cells are refined if the particle mass contained within 
them exceeds a certain threshold value. Thus, this adaptive 
grid follows the collapsing galaxy in a quasi-lagrangian fash- 
ion. With maximum 10 levels of refinement, the peak formal 
spatial resolution is 0.15 comoving kpc. This s imulation is 
described in more detail in lKravtsov et alj j2004). 

Halos and subhalos are identifie d using a variant of the 
Bound Density Maxima algorithm (Klvpin et al. 1999). The 
halo finder algorithm produces halo catalogs, which contain 
positions, velocities, and main properties of the dark matter 
halos between z = 9 and z = at 96 output times, with a typical 



separation between outputs « 10 8 yr. The main properties of 
the halos in the catalogs are: the virial radius (r v ; r ), truncation 
radius (r t ), halo radius (rh a i = rmn[r t ,r v i r ]), halo mass (Mh a i ), 
maximum circular velocity (V max ) and the corresponding ra- 
dius at which the maximum velocity is reached (r max ). The 
virial radius (and corresponding virial mass) is defined as the 
radius within which the density is 180 times the mean den- 
sity of the Universe. The truncation radius is the radius at 
which the logarithmic slope of the density profile becomes 
larger than -0.5, since we do not expect the density profile of 
the halos to be flatter than this slope. If the center of a halo 
does not lie within the boundary of a larger system, its radius 
'"halo is defined equal to the virial radius. On the other hand, 
if the halo is within another halo (and is therefore a subhalo), 
its radius is the truncation radius. The halo mass follows the 
same definition: Mh a i = min[M(r t ),M(r v j r )]. 

In this paper, we consider the evolution of an isolated halo 
Gi, which we call the main halo. This halo contains over 
10 6 particles within its virial radius at z = 0. The virial mass 
at z = in our definition is M v „ = 2.37 x 1O 12 M , the virial 
radius r v ; r = 426 kpc and the halo concentration c = 13. For the 
commonly used overdensity of 340, the virial mass and radius 
are M340 = 2.07 x 1O 12 M and r^o = 330 kpc, respectively. 
The main halo has experienced the last major merger at z ~ 2 
and therefore could subsequently host a large spiral galaxy. 

The main advantage of the catalogs of IKravtsov et alJ 
(2004) is that they trace the progenitors of halos at previous 
epochs and therefore allow us to reconstruct the full dynami- 
cal evolution of the model galaxies. This information is nec- 
essary to follow the evolution of globular clusters within their 
host galaxies, before being accreted onto the main halo. 

2.2. Galactic Potential 

We construct realistic time-dependent galactic potentials 
using the halo catalogs and supplementing each massive dark 
matter halo with a baryonic disk. The distribution of dark 
matter in a cosmological simulation can be approximately de- 
scribed as a combination of isolated halos and satellites within 
them, and therefore the gravitational potential can be approxi- 
mated by a sum of the pote ntials of the main halo and all sub- 
halos ( Kravtso v et al J 120041) . For each galaxy (isolated and 
satellite) we assign a fraction of the mass, Mdisk = /dMmio, 
to the baryonic disk and the remaining mass, (1 -/d)Mh a i , to 
the dark halo. In this crude approximation, the fraction fa ac- 
counts for the potential of stars and cold gas in the galaxy and 
absorbs all the complex details of star formation. We do not 
explicitly consider the galactic bulge because most cluster or- 
bits do not come too close to the bulge, and its potential can 
be described by the inner part of the disk. The bulge mass 
is therefore absorbed in our definition of the disk mass. The 
total potential, as a function of space and time, is the sum of 
the halo and disk contributions of all the galaxies: 

*(r , = ®U(r-rd + <& dlsk (r- r,), (2) 
i 

where r, is the position of the center of mass of halo i at time 
t. 

We assume that the dark matter p otential of each galaxy is 
given by a spherical NFW jNavarro et all 199% profile: 

, , , GM ha i (l ~fd) ln(l +r/r s ) 

$haioM = W1 , r ln , : (3) 

r ln(l + c)-c/(l+c) 

where c = r-^ /r s is the halo concentration, and r s is the scale 
radius approximated as r s m r max /2.16. We use this latter re- 
lation because the profiles of some of the smaller subhalos are 
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not well resolved in the simulation and the radius r max is better 
constrained than r s . The axisymmetric disk potential is mod- 
eled with a Miyamoto & Nagaj] j 19751) profile, in cylindrical 
coordinates (R,Z): 

$disk(^,Z) = j (4) 



in globular clusters in that halo, Mqq, scales with the host 
halo m ass as found in the simulation of iKravtsov & Gnedinl 
< 2005b: 



where Mdisk(0> «d(0 and ba(t) are the mass, radial scale 
length, and vertical scale height of the disk, respectively, 
which vary with time. 

The masses and sizes of the halos are linearly interpolated 
between the output epochs in the halo catalog. For the param- 
eters of the disk, we apply the following simple model. The 
disk mass increases linearly in time from the initial epoch z = 9 
to the redshift of accretion z acc . This is defined as the epoch 
when the halo is accreted onto the main halo, i.e., when its 
center of mass is within the virial radius of the main halo. 
We stop the growth of the disk when it is accreted and keep 
its mass fixed at later times, regardless of the evolution of its 
host halo. For the main halo and other isolated halos, z acc = 
by definition. Thus the initial disk mass is /dMhaioC-Z = 9) and 
the final mass is / d Mh a i (z = Zacc)- We set /d = 0.05, consis- 
tent w ith the value inferred for the Milky Way ( Klvp in et al.l 
2002). 

We use a similar recipe to obtain the scale length «d(/) as 
a function of time, by req uiring q d = 0.02rhain - It is i nterest- 
ing to note that «d for the Miyamoto & Nagaj| (119751) disk is 
equal within 3% to the exponential scale length, r d , defined as 
the radius at which the surface density decreases by a factor 
£(0)/£(>d) = e. The ratio of the vertical scale height to the 
radial scale length is set constant at all ti mes, bdlaA = 0-2, con- 
sisten t with observed galaxy disks (e.g., Binnev & Tremaine 
119871) . In our model, the disk of the main galaxy has the mass 
Mdi s k = 1.2 x 10 n M Q and scale length a d = 8.5 kpc at z = 0. 
Both of these parameters are probably a factor of 2 larger than 
in the Milky Way, as is the total halo mass. 

We use proper (not comoving) coordinates and a reference 
frame where the main halo is at the origin of the coordinate 
system at all times. With the potential given by equation (|2j, 
the acceleration of a test particle at (r,f) is calculated analyti- 
cally: 

a = -V,.$ + ft A Z/ 2 r. (5) 

The second term in this equation is a constant acceleration per 
unit length due to th e dark energy, pa rametrized by the cos- 
mological constant (Laha vet al.l!991l) . It affects only clusters 
outside 200 kpc from the center. We use cubic spline inter- 
polation between the outputs of the halo catalogs to obtain 
the positions of halos at intermediate times. The third order 
spline assures that the halo velocities are smooth functions 
and the accelerations are continuously defined at each catalog 
output. However, this non-linear interpolation of subhalo mo- 
tions creates an additional acceleration of the globular clusters 
within them, which needs to be included in equation (0. 

2.3. Initial Conditions for Globular Clusters 

In the model of Kravtsov & Gnedin (2005), globular clus- 
ters form in dense cores of giant molecular clouds within the 
disks of high-redshift galaxies. We adopt this model to set 
the initial masses, positions, and velocities of the clusters in 
each massive halo, Mh a i > 10 9 M Q . We use the initialglobu- 
lar cluster mass function (GCMF) given by equation ([Q, with 
the slope a = 2 and the normalization such that the total mass 



M cc = 3.2x 10 6 M C . 



alo 



1.13 



x/. 



(6) 



1O U M 0/ 

The additional scaling factor, f > 1, is needed for those mod- 
els where we assume that all clusters form at a single epoch 
rather than over an extended period. We choose this nor- 
malization factor such that the number of surviving clusters 
at z = is approximately equal to the number of observed 
metal poor clusters in the Galaxy. For the best model Sb- 
ii (see below) we have / = 12; the continuous formation 
model Cb-ii does not require a change in the normalization: 
/ = 1 . We adopt the lower limit for the initial GCMF as 
M m m = 10 Mm and the up per limit, M max , given by the fit of 
Kravtsov & Gnedin (2005) scaled by the same factor /. The 
choice of the lower mass limit is motivated by our disruption 
calculations, which show that no cluster of lower mass is ex- 
pected to survive the dynamical evolution over 10 Gyr (see 
§ |4}. The actual number of clusters in each host galaxy is 
given by the discretization of the mass function. 

After being accreted into the main halo, some host halos 
are tidally disrupted. There are a total of 77 halos that host 
globular clusters within the virial radius of the main galaxy, 
of which 37 are disrupted and 40 survive until z = 0. 80% of 
all host halos are accreted in the redshift range z acc w 0.2- 1 .9. 
The surviving halos are accreted systematically later than the 
disrupted halos. The median redshift of accretion for all halos 
is 1.2, while for the surviving subhalos it is 0.64. 

We assume that the spatial number density of globular clus- 
ters follows the distribution of baryon mass in their host 
galaxy disk: dN(R) oc Y, d (R)RdR, where £ d (fl) is the disk 
surface density. We use this equation to obtain the initial ra- 
dial positions of the clusters, R. We assign their azimuthal 
angles and vertical positions randomly from the uniform dis- 
tributions in the intervals [0, 2ir] and [-b d , +b d ], respectively. 
There is no correlation between the initial cluster masses 
and positions. The clusters are set on circular orbits paral- 
lel to the plane of the disk, with the initial velocities equal 
V c = yjGM(R) JR, where M(R) is the total halo plus disk mass 
within radius R. The initial velocities of both halos and glob- 
ular clusters include also the Hubble expansion velocity. 

We consider two sets of cluster formation models: 

1. Single epoch of formation (S). The clusters form at z/ = 
4 in galaxies with Mh a i > 10 9 M Q . The masses and 
positions of clusters are generated using the masses and 
sizes of the halos at z/, but the velocities are obtained 
using the halo properties at the epoch of accretion. In 
the case of the main halo, orbit integration starts at z = 4. 

2. Continuous formation (C). Clusters form continuously 
in the redshift range 9 > Zf > 3, at each of the original 
gasdynamics simulation outputs. The average time in- 
terval between the simulation outputs is 7 x 10 7 yr. The 
span of cosmic time between z = 9 and z = 3 is approxi- 
mately 1.6 Gyr, which is consistent the observed spread 
of ages of the metal-poor clusters in the Galaxy. As in 
model (S), the masses and positions of clusters are ob- 
tained using the halo properties at zj but the velocities 
are calculated only at z acc . 

For each of these formation models, we consider several 
variations of the initial relation between the cluster half-mass 
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radius and its mass and of the evolution of this relation with 
time. For the initial Rh~M relation at the time of formation, 
we consider two models: 

a Rh(Q) = constant. We use the median value for the 
Galactic globular clusters, calculated using the online 
catalog of lHarrisl fi"996): R h = 2.4 pc. 

b Ri,(0) oc M(0) 1//3 . We use the constant average 
half-mass density, Ph(0) = 4 x 10 3 Mp pc~ 3 , from the 
model of Kravtso v & Gnedinl (|2P05) and set R/, = 
(3M/47T/9,,) 1 / 3 . The half- mass density does not depend 
on the position of the cluster in its host galaxy, and 
therefore, by construction does not reflect the local tidal 
field. We make this assumption for consistency with the 
formation model that we are testing. 

For the time dependence of the R/, -M relation, we consider 
three variations: 

i Rh{t ) = constant, 

ii Rh(t) oc M(t )'/ 3 , i.e. p h (t) = constant, 

iii R h (t)ocM(t). 

We construct 12 different model realizations by combining 
the single or continuous formation of clusters (S, C) with the 
initial (a, b) and evolutionary (i, ii, iii) size-mass relations. In 
our notation, Sa-i is the model with a single redshift of forma- 
tion and constant half-mass radius throughout the evolution; 
Cb-i has continuous formation in the redshift range 9 > z > 3, 
initially constant half-mass density, and constant half-mass ra- 
dius thereafter, etc. These models are summarized in Tabled 

For clusters that originate in galaxies that will have become 
satellites of the main halo, we integrate the orbits beginning 
from the epoch of accretion, z acc , rather than the epoch of for- 
mation, Zf- The motivation for doing so is to avoid unnec- 
essary computation of the orbits of clusters inside their host 
galaxies at Zf > z > z acc - Before the host galaxies are accreted 
into the main halo, we expect their globular clusters to remain 
on gravitationally-bound equilibrium orbits. For small hosts, 
the time step required to follow these orbits is small, which 
leads to two problems. One is the large number of compu- 
tations, the other is gradual accumulation of integration er- 
rors. The situation changes when the host galaxy becomes a 
satellite - globular clusters may experience significant tidal 
perturbations from the main halo and may even escape their 
original host. Therefore, we assume that the clusters remain 
bound to their hosts at z > z acc and calculate their orbits only 
from z = Zacc until the present. We keep the cluster positions R 
and Z assigned at zj (equivalent to assuming initially circular 
orbits) but recalculate the initial cluster velocities using the 
host mass at z acc . 

For clusters in all models, we calculate the mass loss via 
two-body relaxation and stellar evolution beginning at their 
assumed epoch of formation, Zf. We include tidal shocks only 
when we begin orbit integration at z acc , because the mass loss 
rate due to tidal shocks is very sensitive to the cluster orbits, 
as described in §|3] 

2.4. Orbits of Host Halos and Globular Clusters 

The orbi ts of clusters in a time-dependent potential pre- 
sented in § 12.21 are integrated using a second-order leapfrog 
scheme with a fixed time step. We have chosen this sym- 
plectic integration scheme because it conserves energy to a 



precision of ~ 10~ 10 in a static spherical potential. A typical 
time step is 10 4 yr, which requires ~ 10 3 integration steps per 
cluster orbit within a satellite host and ~ 10 4-5 steps per orbit 
within the main halo. 

Using only the outputs of an iV-body simulation, we cannot 
reproduce exactly the orbits that test objects would have had 
in that simulation. First, finite integration errors accumulate 
along the orbits and diverge exponentially in time, such that 
two objects initially close to each other may have quite dif- 
ferent final positions. Second, we calculate the gravitational 
potential by approximating the complex distribution of dark 
matter as a sum of spherical halos. This approximation is 
not entirely unr easonable and r eproduces the tidal field fairly 
accurately (see Kravtsov et al. 2004), but of course it is not 
exact. We would like to test the orbit integration scheme on 
objects for which we know their final positions in the simula- 
tion. One class of such objects are the satellite halos that host 
globular clusters. 

We calculate the orbits of the host halos from z acc to the 
present and compare their predicted final radial distribution 
with the actual distribution in the A^-body halo catalogs. In 
addition to the gravitational potential given by equation (|2j, 
we include the effect of dynamical friction, a drag force ex- 
erted by the main ha lo on the orbiting sate llites, using Chan- 
drasekhar's formula dChandrasekharl 19431) : 



ddf = -47rG 2 M ha i p In A - 



erf(X)- 



2X 



(7) 



where Mh a i is the mass of the satellite halo, p(r) is the local 
density of the main halo, v is the velocity of the satellite rel- 
ative to the main halo, and X = v/[2er(r)], where a(r) is the 
one-dimensional velocity dispersion of particles in the main 
halo, which was calculated assuming an isotropi c dispersion 
tensor using the approximate formula of IZentner & Bullockl 
(2003). Assuming a constant value for the Coulomb loga- 
rithm, we find that with In A = 5 the calculated radial distri- 
bution of satellite halos is similar to that in the halo catalogs. 
We somewhat overpredict the number of halos at r < 50 kpc 
and underpredict this number at larger distances, with the 
maximum deviation of ~ 50% at 150 kpc. At even larger 
distances the discrepancy decreases and we match the num- 
ber of halos from the catalogs again at 250 kpc. We find 
that dynamical friction noticeably affects the orbits of satel- 
lit e halos, in an apparent disag reement with the conclusion 
of Penarrubia & Benson (2005). The KS test probability that 
our calculated halo positions are drawn from the same distri- 
bution as the A^-body positions is 23% with dynamical friction 
but only 2% without it. 

Even though we do not reproduce the positions of the sub- 
halos exactly, we find that the calculated orbits are statistically 
consistent with those in the simulation in the region of interest 
(r < 200 kpc). This test therefore suggests that the motions of 
objects in a hierarchically-forming galaxy can be reasonably 
accurately calculated by approximating the total gravitational 
potential as the sum of spherical halos. We now apply this 
method to calculate the orbits of globular clusters. 

Within the disks of progenitor galaxies, all clusters in our 
model begin on nearly circular orbits. Present globular clus- 
ters in the Galaxy could either have formed in the main disk, 
have come from the now-disrupted progenitor galaxies, or 
have remained attached to a satellite galaxy. Figure ^ shows 
the three corresponding types of cluster orbits. Even the clus- 
ters formed within the inner 10 kpc of the main Galactic disk 
do not stay on circular orbits. They are scattered to eccen- 
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FIG. 1 . — Three types of globular clusters orbits. Left panels show the proper (not comoving) distance to the center of the main halo as a function of cosmic 
time, right panels show orbits in the plane of the main disk. Top: cluster formed in the main halo, on an initially circular orbit but was later scattered by accreted 
satellites. Middle: cluster formed in a satellite halo, which survived as a distinct galaxy (halo is shown by a dashed line). Bottom: cluster formed in a satellite 
that was tidally disrupted at t ~ 4 Gyr. 



trie orbits by accreted satellites, while the growth of the disk 
increases the average orbital radius. Triaxiality of the dark 
halo (not included in present calculations) would also scatter 
the cluster orbits. The clusters left over from the disrupted 
progenitor galaxies typically lie at larger distances, between 
20 and 60 kpc, and belong to the inner halo class. Their or- 
bits are inclined with respect to the Galactic disk and are fairly 
isotropic. The clusters still associated with the surviving satel- 
lite galaxies are located in the outer halo, beyond 100 kpc 
from the Galactic center. Note that these clusters may still be 
scattered away from their hosts during close encounters with 
other satellites and consequently appear isolated. 

2.5. Spatial distribution 

Hierarchical mergers of progenitor galaxies ensure that the 
present spatial distribution of the globular cluster system is 
roughly spheroidal. Figure [2] shows the positions of the sur- 
viving globular clusters in the Galactic frame, for the best-fit 
model Sb-ii described in §0] Most clusters are now within 50 
kpc from the center, but some are located as far as 200 kpc. 
The median distance is 90 kpc, which is significantly larger 



than the median distance of the metal-poor Galactic clusters 
(7 kpc). Our simulated galaxy is perhaps twice as extended 
as the Milky Way, but this alone is not enough to explain the 
discrepancy. It could be that some metal-poor clusters have 
formed later within the main disk, which is not included in 
our current model. 

Figure|5]shows that the azimuthally-averaged space density 
of globular clusters is consistent with a power law, n(r) oc r" 7 . 
The best models have logarithmic slopes 7 = 2.65 (Sb-ii) and 
7 = 2.72 (Cb-ii). The slopes for all other models are given in 
Table [0 Since all of the distant clusters originate in progeni- 
tor galaxies and share similar orbits with their hosts (disrupted 
and survived), the distribution of the clusters is similar to the 
overall distribution of dark matter given by the NFW profile. 
The cluster profile is steeper than that of the surviving satellite 
halos (7 w 2.2), because the surviving satellites are preferen- 
tially found in the outer regions of the main halo. A higher 
fraction of the inner satellites is disrupted, but the globular 
clusters that they hosted are dense enough to survive the tidal 
field of the inner galaxy. 
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FIG. 2. — Spatial distribution of surviving clusters in the Galactic frame for model Sb-ii. Dotted circles are at projected distances of 20, 50 and 150 kpc. 
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FIG. 3. — Number density profiles of surviving globular clusters (filled circles) in models Sb-ii (left) and Cb-ii (right) can be fit by the same power law, 
n(r) oc f 2n . Vertical error bars are from Poisson statistics alone. The distribution of model clusters is somewhat steeper than that of surviving satellite halos 
(open circles), but it is similar to the distribution of smooth dark matter (dotted line, arbitrary normalizati on). It is also consistent with the observed distribution 
of metal-poor globular clusters in the Galaxy (solid line), plotted using the data from the catalog of Harri^ 119961) . 



The density profile of model clusters is also similar to the 
observed distribution of the metal-poor ([Fe/H] < -0.8) glob- 
ul ar clusters in th e Galaxy, calculated using the online catalog 
of lHarrisI (1996): 7 = 3.15 ± 0.09. Such comparison is ap- 
propriate, for our model of cluster formation at high redshift 
currently includes only low metallicity clusters ([Fe/H] < -1). 
Thus the formation of globular clusters in progenitor galaxies 
with subsequent merging is consistent with the observed spa- 
tial distribution of the Galactic metal-poor globulars. 

We have also tested the effect of dynamical friction on the 
spatial distribution of clusters in our best models. We included 
this effect on the orbits only within the main halo, but not in 
the original host galaxies. We again use Chandrasekahr's ap- 
proximation for the drag force (eq. 171 ). with the appropriate 
coulomb logarithm, lnA c i = 12, and replace the mass of the 
subhalos by the initial mass of the clusters. By using the ini- 
tial cluster mass, not reduced by subsequent dynamical evolu- 
tion, we consider the maximum effect that dynamical friction 
can have on cluster orbits. Unlike the case of the halo orbits, 



we find that dynamical friction does not noticeably change the 
cluster distribution. The density profile can again be fit by a 
power law with the slope that differs from the case without 
dynamical friction by only A7 = 0.02. This difference is well 
within the errors of the power-law fit procedure and is there- 
fore not significant. We note that the distributions of subhalos 
with and without dynamical friction are also consistent within 
the errors, with A7 = 0.06. 

Several studies in the literature looked at the projected spa- 
tial distribution of globular cluster s in nearby, mainly early- 
type galaxies. Kissl er-Patig l ([1597) find that the surface den- 
sity profiles of globular cluster systems in 15 elliptical galax- 
ies are consistent with a power law a(R) oc R~ r , where 1.0 < 
r < 2.2. More recent studies find similar results for elliptical 
galaxies: T = 1.2 - 1.9 JRhode & Zepfl 12004. T = 0.8-2.6 
( iPuzia et alJ [2004 T » 2.0 jHarris et alJ l2004. T w 1.4 
birsch et all200 5 ),T= 1.2-1.6 ( Forb es et all2006l). and for 
a spiral galaxy NGC 7814, which has T w 1 .0 ( Rhode & Zepfl 
2003). More luminous galaxies have shallower globular clus- 
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FIG. 4. — Mean eccentricity (left) and mean pericenter (right) distributions of the surviving globular clusters in models Sb-ii (top) and Cb-ii (bottom). 



ter pro files, in agreement with an earlier study by lHarrisI 
( 1986). In the inner regions, many GC system profiles have 
constant-density cores. In addition, in all studied cases the 
metal-poor globular cluster population is more spatially ex- 
tended (has lower T) than the metal-rich population. 

To make a prediction for the surface density profile of our 
model clusters, we project the cluster positions along ran- 
dom lines-of-sight and fit a power-law to the resulting distri- 
bution. After averaging over 6000 random projections, we 
find for the two best models: T = 1.76 ±0.11 (Sb-ii) and 
r = 1.81 ±0.12 (Cb-ii). These power-law slopes are consis- 
tent with the above-mentioned wide range of slopes of the 
globular clusters systems in elliptical and some spiral galax- 
ies. These slopes T are also consistent with the expectation 
for a spherical distribution of model clusters, where T = 7— 1 
a nd 7 = 2.7. 

Fort e et all (|2005) emphasize the similarity of the spatial 
distribution of blue globular clusters in giant elli ptica l NGC 
1399 with that inferred for dark matter, and lCote etafl Jl998h 
made a similar case for M87. As Figure [3] shows, the den- 
sity profile of our model clusters is is similar to that of dif- 
fuse dark matter, which is somewhat steeper than the density 
profile of the surviving substructure (e.g., [Nagai & Kravtsov 
12001 iMaccio et alJl2006l iWeinberg et alJl2006l) . A similar 
connection has been made by Moo re et alJ d2006). who identi- 
fied globular cluster hosts with high-density peaks in the dark 
matter distribution. In our model, the orbits of blue clusters 
are determined almost solely by the mass assembly of the host 
galaxy and are therefore only weakly dependent on the final 



morphological type of the galaxy. Therefore, the similar dis- 
tribution of metal-poor clusters and dark matter is expected to 
hold for both spiral and elliptical galaxies. This is a robust 
prediction of our model and is a direct consequence of globu- 
lar cluster formation in progenitor galaxies at high redshift. 

In the continuous formation scenario, clusters that form 
early are more centrally located than those that formed late. 
We split the clusters in model Cb-ii in three equal-size age 
bins and find that the earliest 1/3 of the clusters have a median 
distance of 40 kpc, while the latest 1/3 have a median dis- 
tance of 77 kpc. This result is also a generic prediction of bi- 
ased galaxy formation, where high density peaks that collapse 
early have steeper and more concentrated density profile than 
the m ore common peaks that collapse late (e.g., Moo re et"all 
2006). The oldest metal-poor clusters should therefore be 
found on the average closer to the center of the galaxy than 
the somewhat younger (< 2 Gyr) metal-poor clusters. Note 
that this does not apply to the metal-rich clusters, which prob- 
ably formed later in the main disk of the galaxy. 

2.6. Kinematics 

From the positions and velocities of model clusters we can 
study their kinematics at z = 0. We define orbital eccentricity 
using the pericenter R p and apocenter R a , the radial distances 
from the center of the main halo at the closest and farthest 
points in the orbit, respectively: 



e = 



Kb + Rn 



(8) 
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FIG. 5 . — Anisotropy parameter j3 as a function of radius for models Sb-ii (left) and Cb-ii (right). The vertical error bars represent the error in the mean for each 
radial bin, while horizontal bars show the range of the bin. Horizontal dashed lines illustrate an isotropic (J3 = 0) and a purely radial (J3 = 1) orbital distributions. 



A typical cluster completes several revolutions around the 
center of the main halo, with different R p and R a each time 
because of the time-dependent nature of the potential. We de- 
fine the mean eccentricity, (e), by taking the average over all 
orbital extrema between z acc and z = 0. 

Figure |4] shows the distribution of the average orbital ec- 
centricities and pericenter distances of surviving clusters in 
the two best models. Most orbits have moderate mean eccen- 
tricity, 0.4 < (e) < 0.7, expected for an isotropic distribution. 
The pericenters have a wide distribution, typically between 5 
and 40 kpc from the galactic center. The eccentricities for our 
different models averaged over all clusters are given in Ta- 
ble \l\ They are consistent among different models and lie in 
the narrow range between 0.52 and 0.58. 

The anisotropy parameter is another important measure of 
the overall kinematic properties: 



where (vj:) is the mean square radial velocity of clusters in 
a radial bin, and (v 2 ) is the mean square tangential velocity. 
By construction, (3 = corresponds to an isotropic distribu- 
tion, (3 > is radially anisotropic, and (3 < is tangentially 
anisotropic. 

Figure |5] shows the velocity anisotropy of the model glob- 
ular cluster system as a function of radius in the main halo at 
z = 0, for the two best models Sb-ii and Cb-ii. In the inner 
50 kpc from the Galactic center, (3(r) 0, consistent with an 
isotropic distribution. At larger distances cluster orbits tend 
to be more radial, reaching a maximum value /3 max w 1/2, as 
is typical o f dark matter satellites in iV-body simulations (e.g., 
IColm et all2000l) . There, in the outer halo, host galaxies have 
had only a few passages through the Galaxy or even fall in for 
the first time. The anisotropy of cluster orbits in our model 
derives directly from the orbits of the satellite galaxies. 

The velocity distribution of the model cluster system has a 
moderate rotational component in the inner 50 kpc, V rot w 56 
km s" 1 . The distribution at larger radii is consistent with no 
significant rotation. The plane of rotation coincides with the 
plane of the main disk. Since we do not specify the sense 



of rotation of the Galaxy when we construct a model for the 
gravitational potential, we cannot determine whether the clus- 
ter motion is prograde or retrograde with the respect to the 
disk. 

The alignment of the rotation plane of the inner "disk" 
Galactic clusters with th e plane of the Galactic disk has al- 
ways been expected (Frenk & White 1980; Thomas 1989) but 
could never be confirmed without accurate measurements of 
the cluster proper motions. Our models provide the first in- 
dication that Galactic clusters may indeed rotate along with 
the disk stars. The estimated r otation veloc ity of the metal- 
poor clusters is 30 ± 25 km s -1 ( Harris 200 1 ), which is in rea- 
sonable agreement with V m t = 56 ± 25 km s" 1 for the model 
clusters. 

Figure [S] shows the rotation velocity and the velocity dis- 
persion profile in the model Sb-ii. We calculate the line-of- 
sight dispersion, <ti os , by projecting the 3D velocities along 
10000 random lines of sight and averaging the projected dis- 
persion in six equal-size radial bins. Since the rotation veloc- 
ity component is small, subtracting it has a negligible effect 
on the velocity dispersion, and therefore we do not include 
rotation when calculating the velocity dispersion. The pro- 
jected dispersions, between 80 and 120 km s" 1 , are almost 
identical to the spherical ID velocity dispersions, calculated 
as a\o = <Jt,d I \/3. This is yet another verification that the 
model clusters have an approximately spherical and isotropic 
distribution. 

Model clusters are also fair traces of the gravitational po- 
tential of the dark matter halo. Figure [6] shows two pro- 
jected velocity dispersions for a tracer population with the 
radial density profile of the model clusters, calculated from 
the Jeans equation u sing a publicly available code CONTRA 1 
(Gned in et al.l20 04). The solid line is calculated assuming an 
isotropic velocity distribution, while the dashed line is calcu- 
lated assuming a constant anisotropy parameter (3 = 0.5. If 
the model clusters were in equilibrium with the spherical halo 
potential, in the absence of rotation and of flattening of the 

1 CO NTRA is available for download at 

|http : / / www, astronomy . ohio-state ■ edu/~ognedln/ contra/| 
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FIG. 6. — Velocity dispersion profile of the surviving clusters in model 
Sb-ii. We show the ID (filled squares) and the projected line-of-sight (filled 
circles) velocity dispersions, as well as the rotational velocity of the cluster 
system (filled triangles) in equal-size radial bins. Error bars are obtained us- 
ing bootstrap resampling. Lines show the predicted velocity dispersion in the 
analytical potential of the main dark matter halo, calculated using the spher- 
ical Jeans equation for a tracer population with the same radial distribution 
as model globular clusters, for isotropic orbits (J3 = 0, solid line) and radially 
anisotropic orbits (J3 = 0.5, dashed line). Open circles show the line-of-sight 
velocity dispersion profile of Galactic metal-poor clusters. 

disk potential, we would expect the isotropic distribution to 
reproduce the model ai os at small radii and the anisotropic 
distribution to reproduce it at large radii (see Figure^. Both 
predicted curves agree with the model oi os within the errors, 
but the anisotropic distribution appears to provide a better fit 
at all radii. Figure [6] also shows that the projected velocity 
dispersion of model clusters agrees with that of meta l-poor 
clusters in the Galaxy, calculated using the Harris! i 1996ft cat- 
alog. 

3. DESTRUCTION OF GLOBULAR CLUSTERS 

The mass of a star cluster decreases with time as a result 
of several dynamical processes of internal and external origin 
( Spitzeill987l) . We consider the physical processes that have a 
major effect on the cluster evolution and have been well stud- 
ied in the literature: stellar evol ution , two-body relaxation, 
and tidal shocks. As we note in § !2.5l above. dynamical fric- 
tion does not have a noticeable impact on cluster orbits and 
therefore on their dynamical evolution. 

The differential equation that describes the cluster mass as 
a function of time, M(t), assuming that these effects are inde- 
pendent of each other (a good first-order approximation), has 
the following form: 

dM 

-7- = H>.se + ^ev+fshW, (10) 

at 

where j/ se , v ev , and i/ s h are the time-dependent fractional mass 
loss rates due to stellar evoluti on, two-body relaxa tion, and 
tidal shocks, respectively (e.g.. iFall & Z hang 2001). In the 
following subsections we describe each term of eq. JlOi in 
detail. 

3.1. Stellar Evolution 
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FIG. 7. — Fractional mass loss from a globular cluster due to stellar evolu- 
tion accumulated over time t, for three different IMFs: Salpeter (dotted line), 
Scalo (dashed line), and Kroupa (solid line; adopted in this work). 

We use a simple model to include the ef fect of mass loss du e 
to stellar evolution, following Chernoff & Weinberg ( 1990). 
All stars in the cluster are assumed to form at the same time. 
As stars of different mass evolve away from the main se- 
quence (MS), they lose mass through stellar winds and su- 
pernovae explosions. We assume that the lost stellar material 
is not gravitationally bound to the cluster and escapes imme- 
diately. 

We adopt the stellar initial mass function (IMF) of Kroupa 
d2001l) between O.IM^ and lOOM^: 

dN f mr 13 , m<0.5M Q 
^"Im- 13 , m>O.5M . (U) 

The mass loss in a given range of stellar masses (nt\, m-i) is: 

1*1112 

AM(mi,m2)= / (m — mAdm, (12) 

J mi dm 

where m j is the remnant mass of a star with the initial mass 
m. We adopt t he relation between the i nitial mass and remnant 
mass given bv lChernoff & Weinberg! fl990T) . We use a stellar 
evolution estimate of the MS lifetime for a star of given ini- 
tial mass JChernoff & Weinberg! 1990HHurlev et alJ2000l) and 
convert the mass loss in a range of masses to the mass loss in 
a range of times, &M(m\,m.2) = AM($\,tz). Then we calcu- 
late the fractional mass loss rate of a cluster as a function of 
time, v se (t), with sampling At = 10 7 yr. We have tabulated 
and spline-interpolated these values for use in the dynamical 
calculations. 

Figure shows the fractional mass loss due to stellar evo- 
lution accumulated over tim e, L iA e (t)dt, for alKroup al d2001l) 
IMF and also for a lSalpeterl < 19551) and a lScaloUl986li IMFs, 
which ar e ofte n used in galaxy modeling. For our adopted 
iKroupal J2001I) IMF, a cluster loses « 30% of its mass in 
the first 3 x 10 8 yr because of the evolution of massive stars 
(m > 2M©) and an additional w 10% in the next 10 Gyr due 
to the low mass stars leaving the MS. 

Since the mass loss due to stellar evolution reduces the 
masses of all clusters by the same fraction, this process does 
not affect the shape of the GCMF and only shifts it towards 
lower masses. 
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3.2. Two-body Relaxation 

Internal two-body relaxation, a cumulative long-term effect 
of gravitational interactions between stars that leads them to 
a Maxwellian distribution of velocities, causes some stars to 
gain enough energy to escape from the cluster. The fractional 
mass loss due to two-bod y relaxation ca n be approximated by 
the following expression (Spitzer 1987): 



components: 



£ e 7.25& mG 1 / 2 In A c 



M l / 2 R 



3/2 



(13) 



where £ e is the fraction of stars that escape per half-mass re- 
laxation time, f r h, fh is the mean stellar mass, and In A c i is the 
usual Coulomb logarithm. 

As a good approximation, we adopt constant values for var- 
ious quantities in equation (1131 : £ e = 0.033, the mean value 
obtained fro m A-body and Fokker-Planck models of clus- 
ter evolution (iGnedin et alJU999bl). which is lower than the 
canonical = 0.045 obtained bv lHenonl i 1 96 ll) : m w 0.87 M Q 
for a Rroupa IMF (eq. II II ): lnA c i = 12, a typical value for 
globular clusters, although it is expected to vary somewhat 
with time as the clusters lose mass. With these approxima- 
tions, the time dependence of v ew is in the mass and half-mass 
radius of the cluster: v ev (t) oc M(t)~ x l 2 R h (t)~^l 2 . Since we as- 
sume that Ri,(t) depends only on M(t) but does not depend 
explicitly on cluster's position in the galaxy, the rate of two- 
body relaxation is also independent of the location. This has 
important consequences for the lack of significant radial gra- 
dient of the cluster mass function. 

3.3. Tidal Shocks 

Each time a cluster passes through the disk of its host 
galaxy or near its nucleus, it experiences a rapid change of 
the external tidal force. This change increases the average ki- 
netic energy of stars, reducing the cluster binding energy. If 
the variation of the tidal field (e.g. passage through the galac- 
tic disk) is fast compared with the internal orbital period of 
stars in the cluster, the adiabatic invariants of the orbit are not 
conserved causing some stars to become unbound and escape 
from the cluster. This is more effective in the outer parts of 
the cluster where stars are weakly bound. 

The effects of tidal shocks on the evolution of globular 
clusters have been thoroughly studied in the literature 
( Ostri ker et al]ll972t ISpitzer)ll987t iKundic & Ostrikerll 19951 
I Gnedin & Ostrikerl jl997t Dvlurali & Weinberg! Il997dblal: 
Gnedi net alJ Il999alb!) . The physics is well understood 
and a general mathematical framework has been developed 
to implement shocks due to the flattened (disk-like) and 
spheroidal (bulge-like) components of the galactic potential 
in A-body and Fokker-Planck simulations. These models 
include adiabatic corrections that ensure the conservation 
of orbital invariants when the duration of the shock is long 
compared with the period of stars in the cluster. In such a 
case tidal shocking does not enhance the mass loss. 

We implem ent tidal shocks using the semianalytical theory 
developed by Gned in et alJ (11999a!) and IGnedin! ( 120031) . The 
time variation of the tidal force around a cluster can be di- 
vided into peaks (by absolute value) surrounded by minima, 
and each peak is considered a separate tidal shock. The total 
tidal heating is the sum over all peaks. The tidal tensor is a 
symmetric matrix of the second spatial derivatives of the time- 
dependent potential (eq. 0), including both the halo and disk 



F a =- 



d 2 $ 
dr a drp 



(14) 



A constant level of F a p that does not vary with time is not 
included in the integration; it is considered a steady tidal field 
responsible for tidal truncation of clusters. The ensemble- 
average change in energy per unit mass at the half-mass radius 
can be expressed as 



(A£*) = i(A 



1 



2 6 
where is the tidal heating parameter: 

2 / 



ImR 



hi 



'tid : 



E 

a,0 



-3/2 



F a(} dt ! + ■ 



(15) 



(16) 



'dyn 



where the sum extends over all components of the tidal tensor 
a, f3 = {x,y,z}, in Cartesian coordinates. Here r is the du- 
ration of each tidal shock for each pair of coordinates a, (3, 
and fd y ri is the dynamical time at the half-mass radius. We 
adopt a power-law adiabatic cor rection, appropriate for re la- 
tively extended tidal shocks (see IGnedin & Ostrikerl 19 99 for 
more detail). 

A typical timescale for tidal heating to change the energy 
of stars at the half-mass radius by the order of itself is 

| _ | 
sh (dE h /dt) sh 2{AE h ) ' 
where P is the time interval between tidal shocks, which is 
typically the orbital period of the cluster in its host galaxy. 
The factor of 2 in equation ( 1171 comes from the contribu- 
tion of the second-order energy dispersion, (Aff 2 ), which has 
been shown to be as importan t as the first-order term, (Aff), 
for heating the cluster's stars (Kundic & Ostriker 1995J). The 
energy per unit mass at the half-mass ra dius can be ap proxi- 
mated as \E h \ w vf/2 w 0.2GM/R h (e.g.. lSpitzerll987l) . 

Assuming that a fractional energy c hange in a ti dal shock 
results in the equal fractional mass loss (Chernoff et al. 1986), 
we combine equations (1151 and (1171 to derive the mass loss 
rate due to tidal heating: 



1 



5/3 Aid^A 



(18) 

P GM 

The time dependence of this rate is different from that of two- 
body relaxation: ^ s h(f)ocM(f) _1 ^/,(f) 3 - This difference plays a 
critical role in shaping the mass function of globular clusters. 

Figure [8] shows an example of our calculation of the tidal 
heating parameter, 7 t id- The cluster chosen for this plot was 
formed in a satellite galaxy but was accreted by the main halo 
at t sa 7 Gyr. It is now orbiting the main galaxy and crosses 
the disk every 0.33 Gyr. The dominant component of the tidal 
tensor is F a , which describes the variation of the vertical com- 
ponent of the tidal force during disk crossing. The bottom 
panel of Figure [8] shows the comparison of the tidal heating 
calculated via integration of equation dl6l wit h a traditional 
way o f par ametrizing disk shocks. Following Ostrike ret al.l 
( 119721) and Spitzer ( 1987|, we can write the average energy 
change in one disk crossing as (Aff/,) = 2gl,Rl/3V 2 , where 
g m is the maximum vertical gravitational acceleration dues ot 
the disk, and V z is the component of cluster velocity perpen- 
dicular to the disk. Comparing this traditional expression with 
equation H5\ we find that for disk shocking 

4;e 2 
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FIG. 8. — An example of the calculation of disk shocking. Top panel: 
vertical coordinate of the cluster with respect to the disk of the main galaxy. 
Middle panel: main component of the tidal tensor, F zz , responsible for disk 
shocking. Bottom panel: comparison of the calcul atio n of the tidal heating 
parameter, / t y, using the full tidal tensor (via eq. \16L circles) and using a 
traditional parametrization of disk shocking (via eg. 1191 . triangles). 

As Figure[8]shows, for strong shocks the agreement of the two 
methods is usually very good, whereas for a weak shock the 
traditional expression overestimates the energy change by up 
to a factor of 2. Note that the adiabatic corrections, which we 
omitted here in both cases for clarity, can significantly modify 
the actual amount of heating and should always be included 
in the calculations. 

4. EVOLUTION OF THE GLOBULAR CLUSTER MASS 
FUNCTION 

Using realistic cluster orbits and the specified R/, - M re- 
lations, we now integrate equation dlOi for the evolution of 
GCMF from z = Zf to z = 0. Figure |9] shows the transforma- 
tion of the cluster mass function from an initial power law, 
dN/dM oc M~ 2 , into a final peaked distribution for models 
Sb-ii and Cb-ii. High-mass clusters withM > 2 x 10 5 M Q ap- 
proximately preserve the power-law shape of the IMF, while 
low-mass clusters are more easily affected by the mass loss 
due to stellar evolution, two-body relaxation, and tidal shocks. 

We compare these model mass functions with the distribu- 
tion of metal-poor (MPC, [Fe/H] < -0.8) globular clusters in 
the Galaxy beca use th ey are the type of clusters described by 
our model (see § I2.5> . We obtain V-band luminosities of the 
Galactic clusters from the current version of Harris's catalog 
( Harri sll996l) and use a constant mass-to-light ratio, M/Ly = 3 
in solar units, to transform luminosities to stellar masses. The 
two plotted model mass functions are in excellent agreement 
with the observed GCMF of the Galactic metal-poor clusters. 

We fit the model GCMF with a traditional lognormal func- 
tion, 

dN/dM oc exp[-(logM-logM peak ) 2 /2cr 2 ]rfM. The peak 



mass and the dispersion are very similar for both plotted mod- 
els: M peak = (2.9 ±0.13) x 1O 5 M , a = 0.61 ±0.10 (model 
Sb-ii) and M peak = (2.3 ± 0.13) x 1O 5 M , a = 0.66 ± 0.10 
(model Cb-ii). For the metal-poor Galactic clusters we find 
Wpeak = 2.1 x 10 5 M Q and a = 0.51. The observed and model 
distributions are statistically consistent with each other. 

To further quantify this agreement, we use a Kolmogorov- 
Smirnov (KS) test to calculate the probability of drawing the 
model GCMF and the observed GCMF from the same under- 
lying distribution. Table ^lists the results for all models we 
consider. The two models that we plot in Figure [9] and use to 
study the spatial and kinematic distributions have the highest 
KS probabilities: Sb-ii (P KS = 0.24) and Cb-ii (P KS = 0.06). 
In these models the initial half-mass density is the same for 
all massive clusters and it remains constant through time: 
Rh(t) oc M(Z) 1 ' 3 . A single epoch of formation is slightly pre- 
ferred to the continuous formation scenario, but this differ- 
ence is not statistically significant. In contrast, most of the 
other models are strongly inconsistent with the observed mass 
function of Galactic metal-poor clusters. 

Recall that in setting the initial distribution of model clus- 
ters, we included only massive clusters with M > 10 5 M Q , as 
we expected all lower-mass clusters to be disrupted. We have 
verified this assumption and find that the lowest initial mass 
of the clusters that survive to z = are 3.8 x 10 5 M Q in model 
Sb-ii and 4.4 x 10 5 M Q in model Cb-ii. All of the low-mass 
clusters are destroyed by the present time. 

Table[2shows also the fractions of surviving model clusters 
by mass and by number, for clusters with M > 1O 5 M0. In 
a single formation epoch model Sb-ii, only fy = 0.16 of the 
initial massive clusters survive until z = 0, but they contain 
a higher fraction of the initial cluster mass, /m = 0.46. In 
the continuous formation model Cb-ii, both fractions are a 
factor of 3-4 lower: f N = 0.04 and f M = 0.16. In the latter 
model, low-mass clusters are relatively more abundant at Zf > 
4, when the host galaxies are smaller than at z/ = 4, and their 
subsequent rapid disruption results in the smaller surviving 
fractions. 

Figure [H)] shows the evolution of the cluster mass function 
with time, for the best models Sb-ii and Cb-ii. In model Sb-ii 
all clusters form at redshift z = 4, but in model Cb-ii clusters 
form continuously beginning at z = 9. To illustrate both for- 
mation scenarios in the same way, we plot the distribution of 
clusters in model Cb-ii at redshift z = 4 (formed between z = 9 
and z = 4) as the "initial" GCMF. This epoch corresponds to a 
cosmic time of 1.5 Gyr, which is almost exactly 12 Gyr ago, 
in our adopted cosmology. The sequence of histograms dis- 
plays the progression of the GCMF at four successive epochs 
(0, 1 .5, 6, and 12 Gyr) after z = 4. 

For comparison, Figure [TU| s hows also the m ass functions 
predict ed by the mod el of Fall & Zhang (2001) at the same 
epochs. Fall & Zhane (2001) consider several initial forms of 
the GCMF, several kinematic distributions of the clusters in 
the Galaxy, and several prescriptions for disk shocking, and 
solve equation (II 01 under these different assumptions. They 
also keep the half-mass density M /R\ constant as a function 
of time and find that this is the best way to reproduce the ob- 
served GCMF, in agreement with our conclusions. We select 
one of their models (top left panel of their Fig. 3), which is 
the closest match to our models. This model has an initial 
power-law GCMF with the slope a = 2 and a single epoch of 
cluster formation, which we assume equal to z = 4. We nor- 
malize their model to have the same initial number of clusters 





M (M ) M (M ) 

FIG. 10. — Evolution of the mass function with time, beginning at z = 4 (dotted histogram) and subsequently after 1.5 Gyr (light grey), 3 Gyr, 6 Gyr (dark 
grey), and 12 Gyr (shaded histogram, corresponds to the present epoch). Left p anel: model Sb-ii, ri ght panel: model Cb-ii. For comparison, solid lines (same in 
both panels) show the predicted evolution of the mass function in the model of Fall & Zhans 1 2001) at 1.5 Gyr, 3 Gyr, 6 Gyr, and 12 Gyr after cluster formation, 
from top to bottom. 



with M > 10 5 Mq as in our model in each of the panels of 
Figure ITOl 

We find a reasonable agreem ent of our GCMF with the pre- 
diction of lFall & ZhangH200il) at late times. At 1.5 Gyr after 
cluster formation, our clusters are less affected by the dynam- 
ical evolution, especially in the continuous formation model 
Cb-ii. Note that although our models do not include clusters 
withM < 10 5 M Q initially, this does not affect the comparison 
in the range M = 10 5 -4x 1O 5 M0, where our models predict 
somew hat slower evolution than the model of Fall & Zhang 
teOOll) . The final distributions in our two be st models are 
similar to each other and to the prediction of Fall & Zhang 
J2001I) . The peak of the log-normal fit is logM pea k = 5.46 
(Sb-ii) and 5.37 (Cb-ii) vs. 5.42 (FZ), and the dispersion is 
0.61 (Sb-ii) and 0.66 (Cb-ii) vs. 0.63 (FZ). Note that although 



we use the parameters of the log-normal function for the ease 
of comparison between different models and between models 
and observations, the mass function below the peak is better 
described b y a power law , dN / d log MocMor dN / dM = const 
(see|Fall & Zha ngl2001l) . 

Fi gurefTTIprovides insight into the effects of each of the dis- 
ruption mechanisms separately. Mass loss due to stellar evo- 
lution simply shifts the whole mass function to the left. The 
relative importance of two-body relaxation and tidal shocks 
depends on the Ri,(M) relation. Let us write this relation in 
a general power-law form, R/, oc M s . Then the ratio of the 

mass loss rates is v e v/v s h oc M l / 2 R t ^ 2 oc M 1 / 2 "" 95 / 2 . Unless 
S < 1 /9 (in the plotted model Sb-ii, <5 = 1/3), two-body relax- 
ation is more effective in low-mass clusters, and tidal shocks 
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FIG. 1 1. — Evolution of the mass function in model Sb-ii as a result of 
inclusion of different disruption processes: stellar evolution (light-grey his- 
togram), two-body relaxation and evaporation (dark histogram), and tidal 
shocks (hatched histogram, final distribution). 



are more effective in high-mass clusters. The evolution of the 
Rh(M) relation with time may also affect the balance between 
the two processes. The dark histogram in Figure ITTI shows 
how the inclusion of the evaporation due to two-body relax- 
ation removes a large number of low-mass clusters and turns 
the mass function into a peaked distribution. Inclusion of tidal 
shocks affects mainly the most massive clusters and reduces 
their surviving fraction. 

Our results agree with previous studies of the evolution of 
the cluster mass function, which found that almost any initial 
function can be turned into a peaked distribution as a con- 
sequence of the dynamical evolution. The new result of our 
study is that not all initial conditions and not all evolutionary 
scenarios are consistent with the observed mass function. 

Figure El provides three examples. In the first (top left 
panel, model Sa-i), the half-mass radius is kept fixed at the 
median value for Galactic globulars, R/, = 2.4 pc, for clusters 
of all masses and at all times. The median density M(t)/R\ de- 
creases as the clusters lose mass. Therefore, two-body relax- 
ation becomes less efficient and spares many low-mass clus- 
ters, while tidal shocks become more efficient and disrupt 
most high-mass clusters. The final distribution is severely 
skewed towards small clusters. In the second example (top 
right panel, model Sb-iii), the median density is initially fixed, 
as in our best models, but the half-mass radius is assumed 
to evolve in proportion to the mass, Rh(t) oc M(t). In this 
case the cluster density increases with time. As a result, all 
of the low-mass clusters are disrupted by the enhanced two- 
body relaxation, while the high-mass clusters are unaffected 
by the weakened tidal shocks. The final distribution is skewed 
towards massive clusters. In the third example (bottom left 
panel, model Ca-iii), the initial half-mass radius is the same 
for all clusters but then varies as Ri,{t) oc M(t). In this model, 
which assumes a continuous formation scenario from Zf = 9 
to Zf = 3, the median cluster density again increases with time 
and all of the low-mass clusters are destroyed by two-body 
relaxation. The initial relation Ri t (Q) vs. M(0) does not appear 
to be as important for the final mass function as the evolution 



of this relation with time, Ri,(t) vs. M(t). 

5. RADIAL VARIATION OF THE GLOBULAR CLUSTER 
MASS FUNCTION 

A generic prediction of the dynamical disruption mod- 
els is that the evolution is faster in the inner parts of the 
galaxy, where the stronger tidal field enhances both two- 
body relaxation and tidal shocks. As a result, the low-mass 
part of the mass function should be depleted more strongly 
in the inner regions, causing an appar ent peak mass to b e 
higher than in the outer regions ( Ostri ke7&GnedinllT997h . 
Such radial depende nce of the ma ss function is marginally 
present in the Galaxy ( Gnedin 1997) and more clearly in M3 1 
(Bar mbv et all2001l) . but it has not been detected in other ex- 
tragalac tic systems. This lack of the observed radial variation 
has led Fall & Zhangl J2001I) to consider anisotropic cluster 
orbits, which may suppress radial dependence of the disrup- 
tion processes. 

Our best model Sb-ii does not predict a significant radial 
variation of the mass. To investigate this, we have divided 
the model sample into two subsamples, splitting them at a 
distance r = 30 kpc from the center. The inner sample has 
30 clusters, the outer 70 clusters. (Had we split the model 
sample into two equal parts, the dividing distance would be 
90 kpc, too large to expect any differences in the dynamical 
evolution.) A lognormal fit to the distribution of the inner 
clusters gives the peak mass and dispersion logM pea k = 5.56 ± 
0.20, a = 0.44 ± 0.17, while for the outer clusters logM peak = 
5.39 ±0.17, cr = 0.82 ±0.13, respectively. The KS probability 
of the two subsamples being drawn from the same distribution 
is extremely high, P^s = 0.99. 

The mass functions of the model subsamples are consistent, 
within the la uncertainties, with the mass functions of the 
metal-poor Galactic clusters, split in two equal-size parts at 7 
kpc: logM peak = 5.35 ± 0.09, a = 0.40 ± 0.07 (inner clusters), 
and logMp ea k = 5.31 ±0.16, a = 0.61 ±0.11 (outer clusters). 
Note however the caveat that the assumed dividing radius for 
our model clusters (30 kpc) is significantly larger than that for 
the Galactic clusters (7 kpc). This is partially, but not entirely, 
offset by the fact that our model galaxy is twice the size of 
the Milky Way and the extent of its globular cluster system is 
expected to scale accordingly. 

The lack of a significant radial gradient in our best models 
is due to the contribution of tidal shocks being sub-dominant 
compared to the contribution of two-body relaxation. Under 
our assumption that the median density M/R\ is constant and 
does not depend on the position in the galaxy, the rate of clus- 
ter mass loss due to two-body relaxation alone is independent 
of location. The rate of tidal shocks does depend on the posi- 
tion and on cluster orbits, but because their contribution is rel- 
atively small, the resulting radial gradient is small as well. We 
have verified this hypothesis by looking at one of the failed 
models (Sa-i), where tidal shocks are strongly enhanced com- 
pared to two-body relaxation. This model shows a clear and 
strong radial dependence of the mass function. Therefore, the 
lack of the gradient in our best models may be more of a prod- 
uct of the initial set-up rather than a model prediction, but it 
does show that our best models successfully pass this obser- 
vational test. 

6. CONCLUSIONS 

We have verified the hypothesis that the properties of glob- 
ular clusters formed within high-redshift disk galaxies are 
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FIG. 12. — Models that fail to reproduce the observed mass function of metal-poor globular clusters: with R/,(t) = const (Sa-i, left), and with R),(t) oc M(t) with 
single epoch of formation (Sb-iii, middle) and continuous formation (Ca-iii, right). Normalization of the initial model distribution is such that the number of 
clusters at z = is close to the observed number of clusters in the Galaxy. 



consistent with the observed properties of Galactic metal- 
poor clusters. After calculating the orbits of model clusters 
in the time-variable gravitational potential of a Milky Way- 
sized galaxy, we find that at present the orbits are isotropic in 
the inner 50 kpc and preferentially radial at larger distances. 
All clusters located outside 10 kpc from the center formed in 
satellite galaxies, some of which are now tidally disrupted. 
The spatial distribution of model clusters is spheroidal, with 
a spherical density profile that can be fit by a power law with 
the slope 7 ps 2.7, consistent with the observations of Galactic 
metal-poor clusters. 

Two-body relaxation dominates over other disruption pro- 
cesses and drives the evolution of the cluster mass function 
from an initial power law to the observed peaked distribution. 
The evolution of the mass function in our best models is in 
close agreement with the models of Fall & Zhangl d2001l) at 
late times and is somewhat slower at early times. We find no 
significant variation of the mass function with radius, which 
is largely due to the rate of two-body relaxation being inde- 
pendent of cluster location. Knowledge of realistic orbits is 
important for accurate calculations of the globular cluster evo- 
lution, although we find that tidal shocks play only a moderate 
role in cluster disruption. Therefore, any inaccuracies in orbit 
calculations are not expected to significantly affect the evolu- 
tion of the mass function. 

Interestingly, we find that not all initial conditions and not 
all evolution scenarios are consistent with the observed mass 
function. If the half-mass radius remains constant throughout 
the evolution instead of decreasing with mass, the rate of two- 
body relaxation is suppressed and the final mass distribution 
is skewed towards too many low-mass clusters. On the other 
hand, if the half-mass radius decreases too quickly, propor- 
tionally to the mass, then two-body relaxation is too efficient 
and there are no low-mass clusters left. The successful models 
(Sb-ii and Cb-ii) require the average cluster density, M/R\, to 
be constant initially for clusters of all mass, as in the model of 



Kravtsov & Gnedin (2005), and to remain constant with time. 

We have investigated two formation scenarios. Syn- 
chronous formation of all clusters at a single epoch (z = 4) and 
continuous formation over a span of 1.6 Gyr (between z = 9 
and z = 3) both appear consistent with the observed mass func- 
tion, spatial and kinematic distributions of the Galactic glob- 
ular clusters. However, clusters that formed at later epochs 
have more extended spatial distribution than the clusters that 
formed at early epochs. This is expected for objects that trace 
the distribution of high-density peaks in hierarchical galaxy 
formation. If additional clusters form after z = 3, they are ex- 
pected to have a shallower density profile than the clusters we 
study in this paper. The profile of our clusters is consistent 
with, but is already somewhat shallower than the profile of 
the Galactic metal-poor clusters. Therefore, we do not expect 
a significant fraction of metal-poor clusters to have formed 
after z = 3. Metal -rich clusters, associated with the Galactic 
disk, are not included in our present formation model and will 
be studied in future work. 

In online Tables |2] and [3] we provide catalogs of the phys- 
ical properties of model clusters that survive to the present 
time, for the best models Sb-ii and Cb-ii, respectively. These 
catalogs can be used to compare with other models of the dy- 
namical evolution, to study selection effects in extragalactic 
systems, and to model the kinematics of globular cluster sys- 
tems. 
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TABLE 1 



Model 


R h (0) 


Hh(t) 


7 


w 


r c 
JM 


f d 
JN 


lOgMpeak 


f 

(7 




Sa-i 


const 


const 


2.6 
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0.29 
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1.02 
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Sa-ii 
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0.68 


2 x ur 4 


Sa-iii 


const 


M(t) 
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0.05 


5.98 


0.66 


7 x ur 8 


Sb-i 


M(0)'/ 3 


const 
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4.97 


0.95 


0.0035 


Sb-ii 
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M(0 1/3 
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0.53 
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0.61 


0.24 


Sb-iii 
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0.57 


0.54 
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5.88 


0.52 


4 x 10" 10 


Ca-i 


const 


const 
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0.15 


0.43 


1.96 


1.40 


< io- 10 


Ca-ii 


const 


M(r) 1/3 


2.7 


0.53 


0.11 


0.05 


5.14 


0.71 


0.023 


Ca-iii 


const 


M(t) 
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0.57 


0.16 


0.02 


6.04 
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< io- 10 


Cb-i 


M(O) 1 / 3 


const 
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0.52 


0.17 


0.08 


4.88 


0.83 


3 x nr 4 


Cb-ii 


M(0)'/ 3 


M(0 1/3 


2.7 


0.52 


0.15 


0.04 


5.37 


0.66 


0.063 


Cb-iii 


M(O) 1 / 3 


M(t) 


2.7 


0.56 


0.21 


0.03 


6.08 


0.51 


< io- 10 



"Power law slope of the number density distribution: n(r) oc 

b Time-averaged eccentricity of the surviving clusters. 
c ' d Fraction of the mass and fraction of the number of surviving globular clusters with respect to 
the initial GCMF, for clusters with M > 10 5 Mq . All clusters with the initial mass M < 10 5 Mq 
are disrupted. 

ef Peak mass and dispersion (in log scale) of a lognormal distribution fitted to the mass function 
of model clusters. 

g Kolmogorov-Smirnov probability of the mass function of model clusters being drawn from 
the same distribution as the Galactic metal-poor globular clusters. 
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TABLE 2 

Catalog of clusters in model Sb-ii. 



Zf 


logM 


logM, 


Rh 


Rh,, 


r 
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l°gM ha i ,,- 


X 


y 


z 


Vx 


Vy 


v z 


4.0 


5.72 


6.10 


2.5 


3.4 


4.2 


0.7 


10.9 


-3.8 


-1.6 


-0.6 


-25.4 


-12.3 


-28.5 


4.0 


5.34 


5.91 


1.9 


2.9 


6.4 


l.i 


10.9 


-0.4 


-2.9 


-5.7 


-28.5 


-144.7 


-21.7 
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3.99 


5.72 


0.7 


2.5 


1.4 


1.2 


10.9 


1.1 
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0.7 


57.5 


4.0 


18.6 


4.0 


5.15 


5.87 


1.6 


2.8 


1.4 


1.5 


10.9 


0.8 


0.1 


1.1 


23.1 


112.9 


26.8 


4.0 


6.18 


6.43 


3.6 


4.3 


0.4 


1.8 


10.9 


-0.4 


-0.1 


-0.0 


-110.2 


-3.7 


-94.9 



NOTE. — Columns are: if - redshift of formation, logM - present cluster mass in Mq, logM, 
- initial mass, Rh - present half-mass radius of the cluster in pc, Ri,j - initial half-mass radius, r - 
present distance to the center of the main galaxy in kpc, r, - initial distance to the center of the main 
galaxy, logMhaio,; - mass of the host galaxy at the time of cluster formation in Mq, x,y,z- present 
coordinates of the cluster with respect to the center of the main galaxy in kpc, v x , v v , v z - present 
velocities of the cluster in km s~ . Complete table is available online. 



TABLE 3 

Catalog of clusters in model Cb-ii. 



Zi 


logM 


logM, 




Ri,.i 


r 


n 


logM hato .,- 


X 


y 


z 


Vx 


Vy 


Vz 


6.2 


4.91 


5.83 


1.3 


2.7 


2.3 


1.4 


10.5 


0.7 


2.1 


0.7 


35.3 


32.5 


37.7 


5.6 


4.96 


5.83 


1.4 


2.7 


3.6 


2.8 


10.5 


0.6 


3.5 


1.2 


123.9 


25.4 


64.7 


5.4 


5.01 


5.83 


1.5 


2.7 


2.6 


0.7 


10.6 


-1.3 


-1.9 


-0.7 


-94.4 


-30.2 


-34.0 


5.2 


5.75 


6.13 


2.6 


3.4 


9.8 


1.7 


10.7 


4.2 


7.7 


2.7 


42.5 


7.9 


84.2 


5.0 


5.68 


6.08 


2.4 


3.3 


17.3 


6.5 


10.7 


3.5 


16.5 


4.8 


79.0 


54.0 


10.1 



Note. — 



Columns are the same as in Table|2| Complete table is available online. 



